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Abstract. Given m positive integers R = (r^), n positive integers C = (cj) such 
that = ^ Cj = N, and mn non-negative weights W = {wij), we consider 

the total weight T = T{R,C;W) of non-negative integer matrices (contingency ta- 
bles) D = (dij) with the row sums r^, column sums Cj, and the weight of D equal 

to Y[ '^'ij^ ■ We present a randomized algorithm of a polynomial in N complexity 
which computes a number T' = T'{R,C;W) such that T' < T < a{R,C)T' where 
a{R,C) = min |]^ri!r^ , Y\cj\Cj "^^"^ /N\. In many cases, InT' provides an 
asymptotically accurate estimate of InT. The idea of the algorithm is to express T 
as the expectation of the permanent of an N x N random matrix with exponentially 
distributed entries and approximate the expectation by the integral T' of an effi- 
ciently computable log-concave function on R™"". Applications to counting integer 
flows in graphs are also discussed. 



1. Introduction and main results 

(1.1) Contingency tables. Let us fix m positive integers ri, . . . and n posi- 
tive integers ci, . . . , c„ such that 

ri + . . . + Tm = Cl + . . . + Cn = N. 

A non- negative m x n integer matrix D = (dij) with the row sums ri,... 
and the column sums ci,... ,Crt is called a contingency table with the margins 
R = (ri, . . . ,rm) and C = (ci, . . . The problem of efficient enumeration of 

contingency tables with prescribed margins has attracted a lot of attention recently, 
see [DG95], [D+97], [CD03], [Mo02], [C+05]. The interest in contingency tables 



1991 Mathematics Subject Classification. 05A16, 68R05, 60C05. 

Key words and phrases, contingency tables, permanent, randomized algorithms, log-concave 
functions. 

This research was partially supported by NSF Grant DMS 0400617. The author is grateful to 
Microsoft (Redmond) for hospitality during his work on this paper. 



1 



Typeset by AmS-T^ 



is motivated by applications to statistics, combinatorics and representation theory, 
cf. [DG95] and [DG04]. 

Let W = (wij) be an m X n matrix of non-negative weights Wij. In this paper, 
we consider the quantity 

T(i?,C;V^^) = ^^«'S^ 

D ij 

where the sum is taken over aU contingency tables D — (dij) with the given margins 
R = (ri, . . . , Tm) and C — (ci, . . . , c^). Thus if Wij = 1 for all the value of 
T{R, C; W) is equal to the number of the contingency tables with the given margins. 
If Wij e {0, 1}, the number T(i?, C; W) counts contingency tables D for which we 
have dij = for all i,j with Wij = (here we agree that 0° = 1). In this case, 
T(i?, C; W) can be interpreted as the number of integer flows in a bipartite graph, 
see [B+04] and [C+05]. We note that counting integer flows in a general graph on 
n vertices can be reduced to counting of integer flows in a bipartite graph on n + n 
vertices and hence to counting weighted n x n contingency tables, see Section 1.5. 

Geometrically, one can view T{R, C; W) as the generating function over all inte- 
ger points in the transportation polytope of m x n non-negative matrices with the 
row sums and column sums Cj, cf. [BP99]. 

We note that if m = n, i? = (1, . . . , 1), and C = (1, . . . , 1) then 

T{R,C- W)^peTW 
is the permanent of the weight matrix W, that is, 

n 

pevW = ^Ylwi^(^i^, 

where the sum is taken over all bijections tt : {1, ... , n} — > {1, • • • , n], cf., for 
example, [Mi78] . A randomized polynomial time approximation algorithm to com- 
pute the permanent of a given non-negative matrix was recently obtained by M. 
Jerrum, A. Sinclair, and E. Vigoda [J-l-04]. 

We show that T(i?, C; W) can be represented as the expected permanent of an 
N X N random matrix with exponentially distributed entries. 

Recall that a random variable 7 is standard exponential if 



P(7 >t) = 




for t > 
for t < 0. 



Our starting point is the following result. 
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(1.2) Theorem. Given positive integers ri, . . . , and ci, . . . , c„ such that 

n + . . . + Tm = Ci + . . . + Cn = N 

and mn real numbers Wij, i = 1, . . . , m and j = 1, . . . ,n, let us construct the 
N X N random matrix A as follows. The matrix A = ^(7) is a function of the 
m X n matrix 7 = (7ij) of independent standard exponential random variables ^ij. 
We represent the set of rows of A as a disjoint union of m subsets Ri, . . . , Rm, 
where |-Ri| = for i = 1, . . . ,m and the set of columns of A as a disjoint union of 
n subsets Ci, . . . ,Cn, where \Cj \ = Cj. Thus A is split into mn blocks Ri x Cj. We 
sample mn independent standard exponential random variables ■-yij , i = 1, . . . ,m 
and J = 1, . . . , n, and fill the entries of the block Ri x Cj of A = A{'y) by the copies 
Wij^ij. Then the total weight T{R, C; W) of the m x n contingency tables with the 
row sums ri, . . . , and column sums ci, . . . , c^, where the table D = (Dij) is 
counted with the weight 

wiD) = l[wt;^, 

is equal to 

E per A 
ri! • • -rmlcil ■■■cJ.' 

We prove Theorem 1.2 in Section 2. 

Although we can compute individual permanents per A via the algorithm of 
[J+04], evaluating the expectation is still a difficult problem. However, the expec- 
tation of an approximate permanent of A can be computed efl&ciently. 

(1.3) An approximation algorithm to compute T{R, C; W). We rely heavily 
on the theory of matrix scaling and its applications to approximating the permanent, 
in particular as described in [Lo71], [Si64], [KK96], [NR99], [L+00], and [GS02], as 
well as on the Markov chain based algorithms for integrating log-concave densities 
[AK91], [F+94], [FK99], and [Ve05]. 

We assume here that the weights Wij are strictly positive, which is not really re- 
strictive since the zero weights can be replaced by sufficiently small positive weights. 

Let A = (aij) he an N x N positive matrix. Then there exist positive numbers 
^1, . . . , ^n; ?7i, . . . ,r]N and a positive doubly stochastic (all row and column sums 
are equal to 1) N x N matrix B = B{A), B = (bij), such that 

aij = bij^^rjj for i,j = l,...,N. 

Moreover, the matrix B = B{A) is unique while the numbers ^1,... ,Cn and 
Vij--- jVn are unique up to a scaling 1 — > $^iT, rjj 1 — > r]jT~^ . This allows us 
to define the function 

N 

(t{A) = \[^ir]i. 



We use the two crucial facts about a: 



• There is an algorithm, which, given a positive N x N matrix A and a number 
e > computes (j{A) within relative error e in time polynomial in lne~^ and N 
[L+00] 

and 

• The function a is log-concave, that is, 

Incr {aiAi + 0:2^2) > cii Inu (Ai) + 0:2 Incr {A2) 

for any positive matrices Ai and A2 and any non-negative ai and 0:2 such that 
ai + a2 = 1. 



Our algorithm is based on replacing per A in Theorem 1.2 by the (scaled) function 
o-{A). Namely, we define 



ri\---rm\ci\---cn\' 



Since both o-{A) and the exponential density on M"^"^ are log-concave, and since 
(t{A) is efficiently computable for any positive A, we can apply results of R. Kannan 
et al. [AK91], [F+94], and [FK99] and of L. Lovasz and S. Vcmpala [Ve05] on 
efficient integration of log-concave functions to show that there is a randomized 
fully polynomial time approximation scheme to compute T'{R, C; W). 

• We present a randomized algorithm, which, for any e > computes T'{R, C; W) 
within relative error e in time polynomial in and N (in the unit cost model). 

We discuss the details of our algorithm in Sections 3 and 4. Namely, in Section 
3 we present the necessary results regarding a {A) while in Section 4 we discuss the 
integration problem. 

Finally, we discuss how well the value of T'(i?, C; W) approximates T(i?, C; 1^). 

(1.4) Theorem. For the number T'{R,C;W) computed by the algorithm of Sec- 
tion 1.3, we have 



T'{R, C; W) < T{R, C; W) < a{R, C)T'{R, C; W), 



where 

a{R, C) 




We prove Theorem 1.4 in Section 5. 
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Let us consider the case of m = n and 



n = . . . = Tn = Ci = . . . = Cn = t. 

Thus we are enumerating weighted magic squares, that is, nxn square matrices with 
row and column sums equal to t. Applying Stirling's formula in Theorem 1.4, we 
achieve the approximation factor of a{R, C) < {nt)~^ /'^ {const ■t)'^/'^, that is, simply- 
exponential in the size n of the matrix and polynomial in the line sum t. Thus, for 
any fixed t, the algorithm of Section 1.3 can be considered as an extension of the 
algorithm of N. Linial, A. Samorodnitsky, and A. Wigderson [L+00] for computing 
the permanent of a positive matrix within a simply exponential factor. On the other 
hand, if n is fixed and t grows, the number of magic squares grows as a polynomial 
in t of degree (n — 1)^, see for example, [St97]. Thus, in this case, the algorithm of 
Section 1.3 allows us to capture the logarithmic order of T{R, C; W). 

Apart from the case of Vi = cj = 1 (computation of the permanent), most of 
the research thus far dealt with the case of Wij — 1, that is, with the non- weighted 
enumeration of contingency tables. M. Dyer, R. Kannan, and J. Mount [D+97] 
showed that if the margins are not too small, = (n^m) and Cj = Cl (m^n) , the 
Monte Carlo based approach allows one to approximate the number of contingency 
tables within a prescribed relative error e > in time polynomial in m, n, and 
e~^. In this case, the number of tables is well approximated by the volume of the 
transportation polytope of m x n non-negative matrices with the row sums and 
the column sums Cj. Subsequently, B. Morris [Mo02] improved the bounds to = 
fl (n^/^mlnm) and Cj = (m^/^n Inn). The approximation we get is much less 
precise, but applies to arbitrary weights W = (wij) and seems to be non-trivial even 
for Wij = 1 and moderate values of r^, cj. For example, if m = n and = Cj = n, 
we approximate the number of tables within a factor of (const ■ n)*^"~^)/^, while the 
exact number of tables is at least e'^^"^ \ In other words, in many non-trivial cases 
we get an asymptotically accurate estimate of In T(i?, C; W). 

Since every log-concave density can be arbitrarily closely approximated by the 
push-forward (projection) of the Lebesgue measure restricted to some higher di- 
mensional convex body, the algorithm of Section 1.3 can be viewed as a volume 
approximation algorithm. In contrast to [D-|-97] and [Mo02], the convex body 
whose volume we approximate is not polyhedral. 

(1.5) Counting integer flows in a graph. Let G — {V, E) be a directed graph 
with the set V of vertices and the set E of edges. Hence every edge e e is 
incident to the head head(v) e F of e and the tail tail(e) e V. We assume that G 
is connected and that it does not contain loops or multiple edges. Suppose further 
that each vertex v has an integer number a(f ), called the excess of f , assigned to 
it, and that 
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A set of non-negative integers x{e) : e & E is called an integer feasible flow in G if 
for every v eV the balance condition holds: 



e: head{e)=v e: tail{e)=v 

If G does not contain directed cycles vi vi — > . . . — > ^ 't'l, the number 
of integer feasible flows is finite, possibly 0. The problem of efficient counting of 
integer feasible flows in a given graph has attracted some attention recently, cf. 
[B+04] and [C+05]. A variation of the problem involves introducing capacities of 
edges (upper bounds on the flows). 

One can express the number of integer feasible flows in a graph with = n 
vertices as the number T(i?, C;VF) of weighted n x n contingency tables, where 
Wij G {0, f } for all To this end, let us construct a bipartite graph with n + n 
vertices as follows. For every vertex v E V, we introduce the left copy vl and 
the right copy Vr. The directed edges tt — > v of G are represented by the edges 
ul ^ of the bipartite graph. We also introduce edges vl — > vji. Finally, let us 
choose a sufficiently large integer z, for example, 

z= a{v) 

v. a{v)>0 

and let us assign the excesses 

0' (vl) = z — a{v) and a (vji) = z. 

With a feasible flow in the original graph G we associate a feasible flow in the 
constructed bipartite graph by letting the flow on the edge ul ^ vr equal to the 
flow on the edge u ^ v and assigning the flow vl ^ vr so as to satisfy the balance 
conditions. This correspondence is a bijection between the integer feasible flows in 
G and the bipartite graph. Hence the number of such flows is equal to the number 
of weighted n x n contingency tables with the rows and columns indexed by the 
vertices v E V, the row margins z — a{v), the column margins z and the matrix 
W = (wij) of weights deflned by Wij = 1 for e £^ and Wij = for {i,j) ^ E. 

2. Proof of Theorem 1.2 

If 7 is a standard exponential random variable then for any integer d > we 
have 

r+oo 



Let us consider the random matrix A — {o,pq) as deflned in Theorem 1.2. We iden- 
tify both the set of rows of A and the set of columns of A with the set {1, . . . , N}. 
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For every permutation tt : {1, . . . , N} — > {1, . . . , N}, let 



N 

(2.1) ~ n '^kTv{k) 

k=l 

be the corresponding term of per^. Thus 

(2.2) EperA^^Et,, 

TT 

where the sum is taken over all permutations n. With every permutation tt we 
associate a contingency table D = D{tx), called the pattern of tt as follows. We let 
D = (dij) where dij is the number of indices k E {!,... , N} such that k E Ri and 
7r(/c) e Cj, so the (A;, 7r(A;))th entry of A lies in the block Ri x Cj of A. 
For the corresponding term of the permanent (2.1), we have 

(2.3) Et^ = Y[wf;^dijl 

where D — (dij) is the pattern of tt. 

Now, let us count how many permutations tt have a given pattern D — (dij). 
Let us represent each subset Ri of rows as a disjoint (ordered) union 

n 

i?j = Rij for i = 1, . . . ,m 

3 = 1 

of (possibly empty) subsets R^j with \Rij\ = dij and each subset Cj of columns as 
a disjoint (ordered) union 

m 

^3 = U for j = 1, . . . , n 

of (possibly empty) subsets Cij with |Cy| = dij. This pair of partitions gives rise 
to exactly Y\ij dij\ permutation tt with the pattern D: we choose tt in such a way 
that if A; e Rij then 7r(/c) e Cy and we note that there are precisely dij\ bijections 

Rij ^ ^ij • 

On the other hand, the number of partitions Ri = IJ^ Rij is 

ri\ 

while the number of partitions Cj = (J^ Cy is 
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Therefore, the number of permutations with the given pattern D = (dij) is 



ri! ■ ■ -r^lci! ■ ■ ■ cj 

The proof now follows by (2.3) and (2.2). □ 

Remark. Let us modify the definition of A as follows: instead of filling the Ri x Cj 
block by the copies of Wijjij, we fill Ri x Cj by the copies of just Wij, so A is 
constructed deterministically. It follows from the proof above that the value of 

per A 



ri! • ■■rm^.cil • • - Cn! 



is equal to the total weight of the contingency tables with the margins ri , . . . , 
and ci, . . . ,Cn provided the weight of the table D = (dij) is 



ij J 

(the Fisher- Yates statistics, cf. [DG95]). 

For another proof of Theorem 1.2 in a particular case of Wij = 1, see [Ba05]. 

3. Matrix scaling 

Here we summarize the matrix scaling results that we need. All the results in 
Theorem 3.1 below can be found in the literature 

We reproduce the approach of L. Gurvits and A. Samorodnitsky [GS02] adapted 
to the case of the permanent (paper [GS02] treats a more general and more com- 
plicated setting of mixed discriminants), which is, in turn, a modification of D. 
London's [Lo71] approach. 

Also, we restrict ourselves to the case of strictly positive matrices to avoid dealing 
with certain combinatorial subtleties. 

(3.1) Theorem. For every positive N x N matrix A = (aij) there exist unique 
positive N-vectors x = x{A), y = y{A), and an N x N positive matrix B = B{A) 

X = i^i,--- ,^n) , y = (?7i, . . . , ryjv) , and B = (bij) 

so that the following holds 
(1) We have 



(2) We have 



(^ij = bij^iVj for i,j = 1,... ,N; 

N 

8 



(3) Matrix B is doubly stochastic, that is, 

N N 

= l for j^l,...,N and ^bij = 1 for z = 1, . . . , iV. 
i=i j=i 

Let us define 

N 

Then a is a log-concave function on the set of positive matrices: 

In a {aiAi + q;2^2) > cki In o" (Ai) + q;2 In cr (^2) 

for any two positive N x N matrices Ai and A2 and any two numbers cti, 02 > 
such that CKi + Q!2 = 1. 

Proof. Let us consider the hyperplane 

AT 



H=\{ri,...,rN): J^^i ^ 



in M-^. With a positive matrix A — (aij), we associate the function /a : M-^ 

N / N \ 
fA{t) = In ^ a^jC^^ , where t = (ti, . . . ,tn) ■ 



4=1 



=1 



Then the restriction of fA{t) on H is strictly convex and, moreover, attains its 
unique minimum t* = (r^ , . . . , rj^), t* = t*{A), on H, see [GS02]. 

Since /a is smooth, t* is also a critical point and the gradient of at the critical 
point is proportional to vector (1, . . . ,1), from which we get 



(3.1.1) Yl 



N , ^. 



EN T* 



= 7 



for some constant 7 and A; = 1, . . . , A'". 
Let 

N 

Ci = Y aijc'^o for i = 1, . . . , iV, 

let 

r]j = e~^i" for j = 1,... ,N, 
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and let us define an N x N matrix B = (bij) by 

bij = ^ for iJ = l,...,N. 

We note that 

N 

since t* lies in the hyperplane H with ti + . . . + tjv = 0. Then, by (3.1.1), we have 

TV 

(3.1.2) J]&ii=7 for j = l,...,N 
On the other hand, 

N 

(3.1.3) J2bij = l for i = l,...,N. 

i=i 

Since S is a square matrix, comparing (3.1.2) and (3.1.3), we infer that 7 = 1 and 
so we estabhshed the existence of a; = (Ci, • • • ,Cn) and y = (771, .. . jTJn) and B 
satisfying (l)-(3). 

To show uniqueness, we note that if a; = (^1, ... , ^n), y = {rji, . . . , tjn), and B 
satisfy (l)-(3), then we must have 

N 

6 = 1 = 1,... ,N 

i=i 

and hence, necessarily, the point t = (ti, . . . , tjv) defined by 

Tj = —Inrij for j = 1, . . . , N 

is a critical point of /yi(t) on H. Since /a(^) is strictly convex on H, there is a 
unique critical point t* = t* {A) . 

Thus function a{A) is well-defined. Moreover, we can write 

JV 

(3.1.4) \na{A) = Vln^i = /a {t*) = inmfA{t). 

i=l 

We observe that for any fixed t, the function g{A) = /a(^) is concave on the set of 
positive matrices A = {aij). 
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Hence for any t & H and any ai,a2 >0 such that cti + 0:2 = 1, we have 



/aiAi+aa^alO > (XlfAi{t) + Q!2/yl2(^) > In CJ (Ai) + Q;2ln(J (A2) . 

Taking the minimum over t & H, we conclude that 

Incr {aiAi + 0:2^2) > cti Incr (^1) + 0:2 Incr (^2) , 
so cr(A) is indeed log-concave. □ 

(3.2) Remark. Another useful property of cr(A) which easily follows from (3.1.4) is 
monotonicity: if A — (aij) and A' = (0^^) are positive matrices such that a[j < aij 
for all i and j then cr(A') < a{A). We also note that cr{A) is positive homogeneous 
of degree N: a{XA) = X^a{A) for aU positive N x N matrices A and aU A > 0. 

(3.3) Computing cr{A). N. Linial, A. Samorodnitsky, and A. Wigderson present 
in [L+00] a deterministic polynomial time algorithm, which, given an. NxN positive 
matrix A and a number e > computes the value of ^{A) within a factor of (1 + e) 
in time polynomial in lne~^ and (in the unit cost model). 

We are interested in computing a {A) where A = A('y) is a random matrix of 
Theorem 1.2. Thus A is positive with probability 1. We observe that we can further 
save on computations as follows. 

Let us consider the mxn matrix (wij'jij), which is also positive with probability 
1. Applying the algorithm of [L+OO], we can scale the matrix to the row sums 
and the column sums Cj. Namely, we can compute (approximately, in polynomial 
time) positive numbers Xi,i = 1, . . . , m, and fij, j = 1, . . . ,n, and an mxn positive 
matrix L = (lij) such that 

"^ijlij — hj/J'i^j for i = 1, . . . ,m and j — 1, . . . ,n 
and such that 

n 

hj — fi for i — 1, . . . ,m and 

i=i 

m 

= Cj for j = l,...,n. 

If we divide every row of A from Ri by fiiVi and divide every column from Cj 
by XjCj, we get the NxN matrix with the entries in the Ri x Cj block equal to 
lij/viCj. It is seen that the obtained matrix is doubly stochastic. Therefore, we 
have 

Hence the scaling of the NxN matrix A reduces to the scaling of the mxn matrix 

{Wij^ij). 
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4. Integrating a (A) 
Here we describe an algorithm for computing 

N\ Ea{A) 



T'{R,C; W) = 



ril---rmlcil---CnV 
cf. Section 1.3. 

(4.1) Notation. We interpret the space R"*" as the space of all m x n matrices 
7 = (7jj). Let M!p" denote the positive orthant ^ij > of R"^" and let 

A = I7 : 7y = 1 and 7y > for i = 1, . . . ,m and j = 1, . . . , n| 

be the standard (open) simplex in R"*". 

For < 5 < 1/mn let us consider the 5-interior of A: 

Aj = I7 : 7jj = 1 and > 5 for i = 1, . . . ,m and j = 1, . . . , n|. 



Geometrically, A^ is an open simplex lying strictly inside A. 

For a T e M, let be the Lebesgue measure on the affine hyperplane 



r 



induced by the Euclidean structure on M"^"^. 
For a matrix 7 G W^"^, let 

P(7) = per ^(7) and let ^(7) = a (^(7)) , 

where ^(7) is the matrix constructed in Theorem 1.2 and a is the function of 
Theorem 3.1. 
Thus we have 

T{R,C;W) = ]- / P(7)exp{-^7,,|d7 



and 



T'{R,C;W) = ^ ] , / 5(7)exp{-^7..}rf7, 

+ IJ 



where ^7 is the Lebesgue measure on IR'^"'. 
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To apply the results of [AK91], [F+94], [FK99] (see also [Ve05]) on efficient 
integration of log-concave functions, we modify the problem to that of integration 
of P{'j) and S{'j) first on A and then on As- 

We use that both functions -P(7) and 5'(7) are positive homogeneous of degree 
N and monotone on Rip": if 7 = (7^^) and 7' = (7^'^) are positive matrices such 
that 

l'ij<lij for aU 

then 

P{y)<P{l) and 5(7') < -5(7), 

cf. Remark 3.2. 

(4.2) Lemma. We have 

Proof. We note that 

W^"" = y rA. 

T>0 

Since du dr — sjmn dj, we get 

/ P{j)exp{-J2^,Adj^^ r^e--( [ P{^) dui^)) dr. 

"I" V 

Since -P(7) is positive homogeneous of degree N, we conclude that 
/ P(7) d^.(7) = r^+-"-i / P(7) dz.(7), 

JrA J A 

from which the proof follows. □ 
The same identity holds for the integrals of 5'(7). 

Next, we approximate the integral over the simplex A by the integral over the 
inner simplex A^. 

(4.3) Lemma. Let d < 1/mn be a non-negative number. Then 
(l-mn5)^+'"""' / P(7) dz/(7) < / ^(7) dz/(7) < / P(7) du{^). 

J A J As J A 

Proof. To prove the lower bound, we observe that the transformation 

'jij I — > 7jj — S for i = 1, . . . ,m and j = 1, . . . ,n 
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maps As inside (1 — 6mn)A. Since P is monotone, we get 



P(7) du{-f) > f P(7) du{^) = (1 - 5mn)^+'""-^ / P{-f) du{^), 

J(l-5mn)A J A 



-5mn)A 

where we used that P is homogeneous of degree N . 
The upper bound is obvious. 

The same inequahties hold for the integrals of ^'(7). 
For an < e < 1, let us choose a positive 



□ 



5 < 



-ln(l -e) 



mn{N + mn — 1) 
P(7) dv{-f) 

P(7) dv{-f) 



for small e > 0. 



mn{N + mn — 1) 
Then the integral 

approximates the integral 



within a factor of (1 — e) and the same holds for the integrals of 5'(7). 

Since ^(7) depends linearly on 7, by the results of Section 3, 5'(7) is a strictly 
positive log-concave function on the set of positive matrices 7 and the value of 5" (7) 
can be computed in polynomial time for any given positive matrix 7. 

Our goal consists of estimating T(i?, C; W) by 



T's{R,C;W) = 



1 



(A^ + mn- 1)! 



d^. 



To compute the integral, we apply the algorithms of [AK91], [F+94], and [FK99]. 
The computational complexity of the algorithms is polynomial in the dimension 
mn — 1 of the integral and the Lipschitz constant of InS' on Aj. Hence it remains 
to estimate the Lipschitz constant of In S. 

(4.4) Lemma. Let 6 < 1/mn be a positive number. Let 7 = (7^) and 7' = (7^^) 
be two matrices such that 



Then 



lij,l'ii>S for all 



N 

In S (7) - In 5 (7') I < — max |7y - 7,' • | 

ij 



Proof. For t — {ti, . . . , tn), let 



N / N \ 

flit) = 5^ap9(7)e^' ) , where ^(7) = (apM) 

p=l \q^l J 
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is the matrix of Theorem 1.2. Letting 



by formula (3.1.4), we can write 

(4.4.1) ln5(7) = min^(t). 

Let 

a = max|7y -7 -J. 

Then 

7y < lij + « < 7ij (1 + for all i, j 

and, similarly, 

7ii < lij + a< -fij (1 + for all 

Since 

^pqil) — "^ijlij provided p & Ri and q E Cj, 

we have 

Opq(7) < Op<z(7') (1 + I) 

and, similarly, 

apq(7') < apqil) (1 + f) • 
Therefore, for allt = (ti, . . . , ttv), we have 

A(t) < fj'it) + iVln (1 + I) < fyit) + ^ 

and, similarly, 

aN 

fyit)<f,it) + ^. 

Applying (4.4.1), we complete the proof. □ 

Summarizing, we conclude that there is a randomized algorithm, which, for any 
given e > computes the value of 

nR,C;W) ^^^-'> 



where A is the matrix of Theorem 1.2, within relative error e in time polynomial 
in and N (in the unit cost model). 
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5. Proof of Theorem 1.4 
Our proof is based on two estimates for the permanent of a non-negative matrix. 

(5.1) The van der Waerden bound. Let B — (bij) be an N x N doubly 
stochastic matrix, that is, a non-negative matrix with all row and column sums 
equal to 1. Then 

This bound constituted van B.L. der Waerden's conjecture proved by G.P. Ego- 
rychev [Eg81] and D.I. Falikman [Fa81], see also Chapter 12 of [LWOl]. 

(5.2) A continuous extension of the Minc-Bregman bound. Let B = (bij) 

he an N X N non-negative matrix. Let 

N 

Si = ^bij for i^l,... ,N 

be the row sums of B. 
Ifbij e {0, 1}, the bound 

N 



perS < JJ(si!)^/'* 



1=1 

was conjectured by H. Mine and proved by L.M. Bregman [Br 73], see also Chapter 
2 of [ASOO]. 

A. Samorodnitsky communicated to the author the following extension of the 
Minc-Bregman bound. Suppose that 

N 

(5.2.1) bij = 1 for all for i = l,... ,N 
and that 

(5.2.2) bij <l/ti for all i = l,... ,N, 
where ti,i = 1, . . . ,N, are positive integers. Then 



N 



(5.2.3) perS<JJ 



.=1 



To deduce (5.2.3), we argue that the maximum of perS on the class of N x N 
non-negative matrices satisfying (5.2.1) and (5.2.2) is attained at a matrix with 
bij e {0,1/ti} for all Indeed, let us choose a particular row index i. Then 
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any non- negative matrix B satisfying (5.2.1)-(5.2.2) can be written as a convex 
combination of two non-negative matrices B' and B" which satisfy (5.2.1)-(5.2.2), 
agree with B in all rows, except possibly the zth row, and, additionally, satisfy 
b\j, b'-j G {0,1/ti}. Since the function perS is linear in every row, we conclude that 
per 5 < max{per(i?'), per(i?")}. Proceeding as above for rows i = 1, . . . , N, we 
may assume that bij e {0, for all so (5.2.3) follows from the Minc-Bregman 
bound. 

A similar bound is obtained by G.W. Soules [So03]. If S is a non-negative matrix 
satisfying (5.2.1)-(5.2.2) where do not have to be integer, then 

perB<l[ Y -■ 

Proof of Theorem 1.4- For a given mxn positive matrix 7 = (7ij), let A = ^(7) be 
the matrix constructed in Theorem 1.2 and let B = B{A) be the matrix constructed 
in Theorem 3.1. We have 

per A = a{A) per B, 

and we estimate perB. 

By the estimate of Section 5.1, we have 

W 

pevA>j^a{A), 

from which 

T{R, C; W) >T'{R,C; W). 

On the other hand, as is discussed in Section 3.3, we can construct B = B{^) as 
follows: first, we construct a positive mxn matrix L = (lij) such that 

Wij^ij = kj/iiXj for aU i,j 

and some positive numbers //i, . . . , Hm and Ai, . . . ,n and such that 

n 

lij = ri for i = 1, . . . ,m 

and 

m 

^lij = Cj for j = l,...,n 

and then fill the Ri x Cj block of B by lij/viCj. 

It follows then that every entry of B in the block Ri of rows does not exceed 
l/rj. Applying the bound of Section 5.2, we get 



m , 

r' 
i=i * 

17 



Similarly, every entry of B in the block Cj of columns does not exceed 1/cj, so we 
get 



r ■ 

i=l 



per S < JJ 

Therefore, 

i=l * i=l i J 

and the proof follows. □ 
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